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Abstract 

We derive a formalism of numerical relativity for higher-dimensional spacetimes and develop 



numerical codes for simulating a wide variety of five-dimensional (5D) spacetimes for the first time. 
First, the Baumgarte-Shapiro-Shibata-Nakamura formalism is extended for arbitrary spacetime 



dimensions D > 4, and then, the so-called cartoon method, which was originally proposed as a 
robust method for simulating axisymmetric 4D spacetimes, is described for 5D spacetimes of several 



> 

types of symmetries. Implementing 5D numerical relativity codes with the cartoon methods, we 
{Sj ■ perform test simulations by evolving a 5D Schwarzschild spacetime and a 5D spacetime composed 



of a gravitational-wave packet of small amplitude. The numerical simulations are stably performed 



o 

\ for a sufficiently long time, as done in the 4D case, and the obtained numerical results agree well 



with the analytic solutions: The numerical solutions are shown to converge at the correct order. 
We also confirm that a longterm accurate evolution of the 5D Schwarzschild spacetime is feasible 
using the so-called puncture approach. In addition, we derive the Landau-Lifshitz pseudotensor 
in arbitrary dimensions, and show that it gives a robust tool for computing the energy flux of 
gravitational waves. The formulations and methods developed in this paper provide a powerful 
tool for studying nonlinear dynamics of higher-dimensional gravity. 

PACS numbers: 04.25.D-, 04.50.-h 
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I. INTRODUCTION 



Clarifying the nature of higher-dimensional gravity has become an important issue, since 
the braneworld scenarios were proposed If the space in which we live is a three- 

dimensional (3D) brane in extra spatial dimensions that are large or warped, the Planck 
energy may be of O(TeV) and quantum gravity phenomena may emerge in high-energy 
particle colliders such as the LHC. If this scenario is correct, mini black holes may be 
produced at the LHC , and this fact motivated a lot of theoretical works in the past 

decade (see |6| for a recent review). Understanding the AdS/CFT correspondence is also an 
interesting issue in the higher-dimensional gravity. 

To study nonlinear dynamics of spacetimes, numerical relativity is probably the unique 
approach. In the past decade, numerical relativity for four-dimensional (4D) spacetimes 
was significantly developed. Now, it is feasible to perform a longterm simulation for merger 

ack holes, which is one of the 



of binary black holes or for high-velocity collision of two b 

strongest gravitational phenomena in nature (see Refs. 0, y, y, y, I 111 . Il2j ] for pioneer works 

of binary black hole merger). It is natural to expect that the formulation and numerical 
techniques developed for 4D cases can be extended to the higher-dimensional cases. 

There are also a few pioneer works in the five- dimensional (5D) numerical relativity 
performed in the past decade Q • However, the purpose of these works was to study a 
specific issue, i.e., the Gregory-Laflamme instability of a black string. Thus, the formulation 
and numerical method in these works are applicable only for this particular issue, and thus, 
developing a general formulation and codes in the higher-dimensional numerical relativity is 
still an issue. Furthermore, there obviously remain a lot of interesting issues to be explored 
in this field, as partially listed in the following. 

The first issue is black hole formation in high-energy particle collisions. If a black hole 
is formed at the LHC, it will emit the Hawking radiation and may be detected. To predict 
the rate of mini black hole production and its detectability, it is necessary to know the cross 
section for the black hole production o"bh and the resulting mass and angular momentum 



of the formed black hole. A partial answer was given in Refs. [la, |16[] (see also [Hj) by 



numerical 
particles 



y 



solving the apparent horizon at an instant of the collision of Aichelburg-Sexl 



18 1 in higher dimensions. Because the apparent horizon formation implies the 

n 

formation of the event horizon assuming the cosmic censorship (e.g., see [19(), the cross 
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section of the apparent horizon formation a ah gives the lower bound of <7bh- However, 
the precise value of <tbh itself is necessary for exactly predicting the phenomena in the 
particle collider. In the 4D case, high- velocity collisions of two relativistic objects have been 
studied by full numerical relativity, via a model of high- velocity collision of two black holes 



20, 



21 



221 ] . In particular, Ref. 



211 ] (see also Ref. [22| for the refinement of the work of 



2 1! ]) first studied the collisions with nonzero im pac t parameters and clarified that o"bh is 



approximately twice as large as a ah found in Ref. [16j. They also showed that resulting mass 
and angular momentum should be significantly smaller than the initial values of the system 
because of a huge amount of gravitational radiation. However, these studies are nothing but 
a prelude of the study of high-velocity collisions in higher-dimensional spacetimes, which is 
really required. 

The second issue is on the stability of higher-dimensional rotating black holes (Myers- 
Perry black hole) 23j. Although there are works on the stability of those black holes by 
separating variables for the metric perturbation in the linearized Einstein equation, the 
anal ysis can be applied for limited situations (see, e.g., 24| for special rotation parameters 
and [25[] for a tensor-mode perturbation) and the problem has not been entirely investigated. 
Hence, multidimensional numerical analyses are required. Clarifying the stability of the 
higher-dimensional rotating black holes with a single rotation parameter, a, is important, 
because such a black hole would be the outcome of particle collisions in the TeV gravity 
scenarios unless the formed black hole is unstable. However, a recent numerical analysis of 
the linearized Einsteins equation suggests that black holes of high values of a are unstable 
at least for spacetime dimensions D > 7 [26!] (see also 2^]): The production rate of mini 
black holes may be smaller than what we naively expect from the analysis of the apparent 
horizon jig . To elucidate the stability and subsequent evolution after the onset of instability, 
numerical relativity will play a crucial role. 

The third issue is on the evolution of a black hole on a Randall-Sundrum (RS) brane. 
So far, no analytic solution of a 5D static black hole localized on the RS brane has been 
found. The recent numerical work 28j| (see also 29() indicates the nonexistence of such 



solutions. If this is the case, any black hole produced on the RS brane cannot relax to a 
stationary state but evolve in time. Clarifying the fate of such black holes is an interesting 
issue. Furthermore, if the AdS/CFT correspondence holds for the RS models, a 5D classical 



black hole on the RS brane is dual to a 4D black hole with quantum fields 



30, 



31] , and thus 
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we could obtain an indication for the Hawking radiation including the back-reaction effects. 

Motivated by these issues, we have developed numerical relativity codes for simulating 
5D spacetimes as the first step. The purposes of this paper are the following three. The first 
purpose is to describe a numerical relativity formulation in higher-dimensional spacetimes. 



We adopt the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism 32|, |33| and write 
down its higher-dimensional version. 



The second purpose is to describe the cartoon method [34], [35| for 5D spacetimes of 
several types of symmetries. The cartoon method was originally proposed for simulating 
4D spacetimes of axial symmetry using the Cartesian coordinates. This method has been 
demonstrated to be quite robust for accurately and stably simulating not only yacuum 
spacetimes but also rotating stars and rotating stellar core collapses (e.g., Ref. |36|). The 
essence of this method is that we do not have to use curvilinear coordinates that have 
coordinate singularities. In most higher-dimensional problems, the spacetime should have 
symmetries, e.g., among the extra-dimensional directions. For such problems, it will be 
better not to adopt the curvilinear coordinates but to adopt the Cartesian coordinates for 
an accurate and stable simulation, as we have learned in the 4D simulations. In the higher- 
dimensional issues that we listed above, several types of symmetries may be imposed. For 
example, in the off-axis collision of two black holes, the axes perpendicular to the orbital 
plane should be equivalent. In this paper, we particularly focus on 4D spaces (i.e., 5D 
spacetimes) with U(l) symmetry, U(l) x U(l) symmetry, and SO(3) symmetry. 

The third purpose is to report our new codes for simulating 5D spacetimes, which are 
implemented using the BSSN formalism and cartoon methods. For demonstrating that the 
codes work well, we perform simulations for test problems for which analytical solutions are 
known. Specifically, evolutions for a 5D Schwarzschild spacetime and for a 5D spacetime 
composed of a gravitational-wave packet of small amplitude are chosen for the tests. In 
the former case, we first solve the 5D Schwarzschild spacetime by our codes in the geodesic 
slice and show that the numerical results agree with the analytic solution derived in this 



paper. We also evolve the spacetime by the puncture approach [8( with the dynamical 
slices and T-driver conditions K 



and show that the longterm evolution of a black hole 
spacetime is feasible as in the 4D case. In the second test, we compare the numerical 
results of a gravitational-wave packet with the semianalytic solution for linearized Einstein 
equations given in Appendix A, and show that they agree well. In addition, we study a 
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method for estimating the energy flux carried by gravitational waves. The Newman-Penrose 
formalism is widely used for extracting gravitational waves in the 4D numerical relativity. 
Unfortunately, such formalism has not yet been developed in hig her dimensions. Here, we 
propose the higher-dimensional Landau-Lifshitz pseudotensor |40j for calculating the energy 
flux carried by gravitational waves and demonstrate that it correctly gives the amount of 
radiated energy. 

This paper is organized as follows. In Sec. II, we derive the BSSN formalism in higher 
dimensions. In Sec. Ill, we describe the cartoon methods in 5D spacetimes of the three types 
of symmetries listed above. In Sec. IV, we present the numerical results of test simulations 
for 5D numerical relativity and show that they agree with the analytic solutions. We also 
show that the energy extraction by the Landau-Lifshitz pseudotensor works well. Section V 
is devoted to a summary. In Appendix A, we summarize the equations and analytic solutions 
for the linearized 5D Einstein equations of U(l) x U(l) symmetry or 5*0(3) symmetry. The 
Landau-Lifshitz pseudotensor in higher dimensions is given in Appendix B. Throughout the 
paper, we use the units of c = 1 where c is the speed of light. 



II. BSSN FORMULATION IN HIGHER DIMENSIONS 



32 



331 ] for higher-dimensional space- 



In this section, we describe the BSSN formalism 
times. After reviewing the Arnowitt-Deser-Misner (ADM) formulation for D dimensions in 
Sec. IIA, the D- dimensional BSSN formalism is derived in Sec. IIB. 



A. ADM formulation 



Suppose M. be a D-dimensional spacetime with a metric g a b- Consider a sequence of 
(D — l)-dimensional spacelike hypersurfaces £4(7^, K ao ) foliated by a time coordinate t in 
M.. Here, 7^ is the induced metric of S t defined by jab '■= 9ab + n a n b with the future- 
directed unit normal n a to E t . K ab is the extrinsic curvature K ab := — (l/2)£ ra 7 afe , where C n 
is the Lie derivative with respect to n a . The coordinate basis t a of the time coordinate t is 
decomposed as t a = an a + /3 a , where a is the lapse function and (3 a is the shift vector. 

The .D- dimensional Einstein equation ( D 'G a t, = SnGoTab is decomposed into constraint 
and evolution equations. Here, ^G a b, T a b, and Gd are the D- dimensional Einstein tensor, 
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the stress-energy tensor, and the gravitational constant, respectively. First, we define 

p := T ab n a n b - 3a := -T 6c nV a ; S ab := T cdl c al d b . (1) 



The Hamiltonian constraint is derived from the Gauss equation to give 

R + K 2 - K ab K ab = lQnG D p, (2) 

where R is the Ricci scalar of the spacelike hypersurface S t . The momentum constraint is 
derived from the Codacci equation to give 

D b K\ - D a K = 8nG Dja , (3) 

where D a denotes the covariant derivative with respect to •y ab . The evolution equation of 
the induced metric 7^ is derived from the definition of the extrinsic curvature as 

Ctlab = -2aK ab + D a f3 b + D b f3 a , (4) 

and the evolution equation of the extrinsic curvature K ab is derived from the Ricci equation 
to give 



C t K ab = -D a D b a + a {R ab - 2K ac K c b + K ab K) 



+ (3 c D c K ab + K cb D a f3 c + K ca D b f3 c - %itG D a 



dab T jy 2^ ab 



(5) 



where R ab denotes the Ricci tensor with respect to j ab and S := S c c . The D-dimensional 
equations are formally different from the 4D equations only in the coefficient of the last term 
of Eq. In vacuum, the equations are independent of the value of D. 

The above expressions are given in the covariant way. Here, we introduce coordinates 
x l that span the hypersurface S f , where i,j = 1, ...,£) — 1. In these coordinates, the line 
element is written by 

ds 2 = -a 2 dt 2 + jijidx* + ftdt^dxi + (3 j dt), (6) 

and the spatial components 7^ of 7°* are the inverse of 7^. The constraint and evolution 
equations in the coordinate expressions are obtained just by replacing the indices a, b to 
the spatial indices i, j and the Lie derivative C t to the coordinate derivative d t in Eqs. (J2J), 
©, ©, and ©. 



B. BSSN formalism 



Now, we derive the BSSN formalism for higher- dimensional spacetimes. The basic idea 
of the BSSN formalism is to increase the number of variables as well as that of constraints in 
order to guarantee the stability in numerical computation (e.g., to kill constraint violation 
modes). First, 7^ is conformally transformed as 

Jij = Xliji (7) 

where the conformal factor x is chosen so that the determinant of 7^ (denoted by 7) satisfies 
the condition 

7 = 1. (8) 

This is equivalent to setting x = 7~ 1 /( D ~ 1 \ We choose 7^ and x as the fundamental 
variables. 

In the original BSSN formalism for the 4D spacetime, the conformal factor e -4 * was used 
rather than x- I* 1 t ne 4D puncture formalism, x or W = x 1 ^ 2 is often used js]. For evolving 
the puncture black holes in five dimensions, x turned out to be a good choice. This is the 
reason that we choose x as one of the fundamental variables. 

Next, the extrinsic curvature is decomposed into the trace part and the trace- free part as 

K 

Kij = Aij + — — 7y, (9) 

where K denotes the trace of and Ay is the trace-free part. As in the 4D BSSN 
formalism, K is chosen to be one of the fundamental variables. The trace-free part is 
conformally transformed as 

Aij := xAij, (10) 

and A^ is chosen to be one of the fundamental variables. Hereafter, the indices of A^ and 
A 13 are raised and lowered by the conformally transformed metric ^ and 7V/- 

In terms of the variables Xi 7^-, and A^, the Hamiltonian constraint ([2]) and the 
momentum constraint are rewritten as 

R + ^j k2 ~ h~ AlJ = levrG^p, (11) 

and 

di& + r{ k A ik - Bjllfi K . _ = S7rG D j ij j h (12) 

L> 1 Z x 



where T J ik is the Christoffel symbol with respect to 7^ and the comma (, i) denotes the 
derivative by x l . 

The evolution equation of \ is derived from Eq. (j3J) with Eqs. (j7j) and (jSJ) to give 

(d t - pfyx = (ocK - d^) . (13) 

Multiplying ■y ab to Eq. flSJ and rewriting it with Eqs. ([2]), (TTOT) . and (ITTj) . the evolution 
equation of .K" is derived to give 

(d t - f3%)K = -DiD'a + a ^ A tl + + [(£> - 3)p + 5] . (14) 

Rewriting Eq. PJ with Eqs. ([7]), (JUJ), (TTOT) . and (TT5|) . the evolution equation of the conformal 
(L> — l)-metric is derived as 

{dt - (3 k d k )% = -2aA, + 7^/3* + % k d^ k - -^jd k (3% r (15) 

The evolution equation of A\j is derived by substituting Eq. with Eqs. (J7j) and (fT0|) into 
Eq. © and using Eqs. (fTT]h (fT3]h (fH|h and ([T5]), to give 

[dt - /M)4 = X [-(A^«) TF + a (i?J F - 8vr5j F )] 

+ a (kI^ - 2A ik A]) + 4*^/3* + i^/J* - -^—^ k ^ k A^ (16) 

where TF denotes the trace-free part, e.g. Rj? = Ry — R^y/(D — 1). 
The Ricci tensor is decomposed into two parts as 

Rij = Rij + (17) 

where R^ is the Ricci tensor with respect to %j and R^ is the contribution of the conformal 
factor. Here, R^ has the terms (1/2)7^(7^ + ju tk j ~lki,ij ~lij,ki) an< ^ ^ ne ^ rs * three terms 
could be the source of a numerical instability, as often found in 4D numerical relativity. For 
stable numerical integration, the following auxiliary variable is introduced 33]: 



y*fj* = -7%- (is) 



We note that another choice := S^ k d k jij can be used as well 32]. It was found that the 
numerical results for the test simulations in this paper do not essentially depend on the 
choice. 
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~ ' ~ (y) 

Using the variable P, Rij and R\f are rewritten as 



R 

and 
R 



ij 



-pl "pA; 
1 ik L jb 



(19) 



(x) _(D-3) 



2X 



X,i, 



P - 3) x> a, 



4 



2 X 



(20) 



2% v _/ 4 X 2 

where 7 = 1 is used in deriving these equations. As in the 4D case, the second derivatives 
of %j explicitly appear only in the first term of Eq. (|19|) . 

Since P is one of the dynamical variables in the BSSN formalism, its evolution equation 
has to be derived. Substituting Eq. (Tl5|) into d t T l = dj(Y k l' il lki,t) and eliminating A 13 ■ 
with Eq. ffl2|) . we obtain 

(d t - pdj)r = -2A ij d ja + 2a 



1 

2 



(D-1) X 



l A ij 



D-3 



2 X 



(21) 



D-l Jl D-l 

In summary, the variables to be evolved are x> 7y, A^- and T* (or Fj), and their 
evolution equations are Eqs. (TT5|) . (1T41) . (TT5|) . (TTBT) . and (|2T|) . respectively. The conditions 
(JH]), tr(Aij) = 0, and (fl8i) are regarded as the new constraints which arise because the 
dynamical variables are increased. 

As shown above, the BSSN formalism for higher dimensions essentially has the same form 
as that for the 4D case, except that some coefficients are changed. Because of the change 
in the coefficients, the behavior of the solutions near the black hole and in the wave zone 
is changed. However, this change does not significantly affect the stability and accuracy in 
numerical computations at least for evolutions of the 5D Schwarzschild spacetime and a 5D 
spacetime of small- amplitude gravitational waves as shown in Sec. IV. 



III. CARTOON METHOD 



In this section, we describe the cartoon method for 5D spacetimes. The cartoon method 
was originally proposed as a prescription for stable numerical simulations of axisymmetric 4D 

9 



spacetimes. The essence in this method is not to use curvilinear coordinates but to use the 



Cartesian coordinates 



341 ] . We briefly review this (say, the case u x = y, z") in Sec. IIIA. 



Next, we extend this method to 5D spacetimes with symmetries. In higher-dimensional 
spacetimes, there are various types of possible symmetries. Here, we consider 4D spaces 
(i.e., 5D spacetimes) with U(l) symmetry (the case u x, y, z = w"), U(l) x U(l) symmetry 
(the case u x = y, z — w"), and 5*0(3) symmetry (the case u x = y = z, w"). The cartoon 
methods for these three cases are described in Sees. IIIB, IIIC, and HID, respectively. 



A. 3D axisymmetric space 

For 3D axisymmetric spacelike hypersurfaces in a 4D spacetime, the 3D Cartesian coor- 
dinates (x,y,z) can be introduced so that the vector <9 V := xd y — yd x becomes the Killing 
vector. In other words, each spacelike hypersurface has U(l) symmetry around the z axis. 
We refer to this case as x = y, z in short. 

One natural coordinate choice for this space is the cylindrical coordinates (p,ip,z). If 
these coordinates are adopted, the problem reduces to a 2D problem (i.e., all quantities 
depend only on p and z). However, in these coordinates, the symmetry axis p = is the 
coordinate singularity. On this coordinate singularity, one has to change the manner of 
finite differencing because there is no point of p < 0. This sometimes (not always) causes 
a numerical instability, which is known as the finite discretization instability. Although it 
might be possible to stabilize numerical computations by appropriately modifying the finite- 
differencing method, there is the case that the prescription is not simple or has not been 
found without numerical viscosity, e.g., issues for which a longterm simulation of rotating 
objects is necessary. 

One can avoid this problem by using the Cartesian coordinates because they have no 
coordinate singularities. The shortcoming in the Cartesian coordinates is that U(l) sym- 
metry does not explicitly appear in equations, and thus, we have to solve 3D equations. 
Suppose that the initial data are given on the (x, z) plane (i.e., <p = 0). In the case that the 
cylindrical coordinates are adopted, the subsequent evolution of the system is feasible with 
this data. However, in the Cartesian coordinates, one cannot calculate the next step only 
with this data, because the equations include y derivatives of functions to be solved. 

However, we do not have to prepare the data for all values of y, if the cartoon method 
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is used. In this method, a few grid points in the neighborhood of the (x, z) plane are 
prepared. The number of necessary grid points depends on the order of numerical accuracy 
required in the finite differencing (see the last paragraph of this subsection). Then, the 
data at a grid point (x,y ^ 0, z) are generated using the data at a point (p, 0, z) [i.e., on 
the (x,z) plane], where p = \J x 2 + y 2 , by use of the symmetry. Because the grid point is 
not located at the point (p, 0, z) in general, the data at this point are determined by an 
interpolation. The method of the interpolation depends on the required order of accuracy 
(see the last paragraph of this subsection). Once the data at the grid points y ^ are 
known, y derivatives at y = are calculated and the data on the (x, z) plane are evolved 
toward the next time step. 

The symmetric relations are derived by the fact that the Lie derivative of functions with 
respect to the Killing vector becomes zero. For a scalar function ty(x,y,z), the symmetric 
relation is 

y(x,y,z) = y(p,0,z). (22) 

In order to derive the symmetric relation of a contravariant vector function T l , it is conve- 
nient to consider the coordinate transformation from the (p, <p) coordinates to the Cartesian 
coordinates. After expressing T x (x,y,z) and T y (x,y,z) in terms of T p (p,z) and T v (p, z), 
the latter two can be replaced by the relations on the (x, z) plane, T x (p, 0, z) = T p (p, z) and 
Ty(p, 0, z) = pT^{p, z). This yields 

T*{x, y, z) = (x/p)T x (p, 0, z) - (y/p)T^(p, 0, z), (23) 
V(x, y, z) = (y/ P )T x (p, 0, z) + (x/p)T^(p, 0, z). (24) 

The relation between T z (x, y, z) and T z (p, 0, z) is the same as that for a scalar, described in 
Eq. fT22l . A covariant vector Ti(x, y, z) has the same symmetric relation as that of T l (x, y, z). 

In a similar manner, the symmetric relation of a symmetric covariant tensor function 
Sij = Suj) is obtained. S zz (x,y, z) has the same relation as that for a scalar, Eq. (T22J), and 
S zx and S zy have the same relations as x and y components of a vector function, Eqs. (I2"3"j) 
and fT24"|) . For the other components, the following relations are derived: 

S xx (x,y,z) = (x/p) 2 S xx (p,0,z) + (y/p) 2 S yy (p,0, z) - (2xy/ p 2 )S xy (p,0, z), (25) 
S yy (x,y,z) = (y/p) 2 S xx (p,0,z) + (x/p) 2 S yy (p,0,z) + (2xy/p 2 )S xy (p,0,z), (26) 
S xy (x,y,z) = (xy/p 2 )[S xx (p,0,z) - S yy (p,0,z)} + [(x 2 -y 2 )/p 2 ]S xy (p,0,z). (27) 
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Again, a contravariant symmetric tensor has the same symmetric relation as that of Sij. 

Using the above relations, the data for y ^ are generated using the data in the (x, z) 
plane, and thus, the derivatives with respect to y can be calculated. The required grid 
number is 5 for the fourth-order finite differencing (i.e., the data at y = ±Ay and ±2 Ay 
have to be determined, where Ay is the grid spacing), and 3 for the second-order one. For 
obtaining the values at a point (p, z) on the (x, z) plane, interpolation is necessary. To keep 
the fourth-order accuracy, we have to use at least fourth-order accurate interpolation (e.g., 
fourth-order Lagrangian interpolation). 

B. 4D space with U(l) symmetry 

In the following, we describe three cartoon methods in 5D spacetimes (4D spaces) of 
three types of symmetries, denoting the Cartesian coordinates by (x, y, z, w) and assuming 
that the 4D space is topologically identical to the 4D Euclidean space. 

First, we consider a 4D space of U(l) symmetry whose corresponding Killing vector is 

= zd w — wd z (i.e., tan?/; = w/z). An example of a system of this symmetry is an off-axis 
collision of two black holes. Suppose that the centers of the two black holes are located in 
the (x, y) plane. In this case, the directions orthogonal to the (x, y) plane (i.e., the direction 
of z and w axes) are equivalent, and thus, the space has the U(l) symmetry. We refer to 
this symmetric space as x, y, z = w in short. 

Such a spacetime can be simulated as a 3+1 problem using the cartoon method in a 
similar prescription to that in the 3D axisymmetric space. We first prepare the grid points 
in the (x, y, z) plane and a few neighboring grid points in the w direction. Then, the data 
at a point (x, y,z,w ^ 0) are generated by the data at a point (x, y, p, 0) with symmetric 
relations, where p = \J z 2 + w 2 . The symmetric relations are essentially same as those in 
the 3D axisymmetric case: It is sufficient to replace the indices (x, y) in Eqs. (12"3"1) - (|2"TI) to 
(w, z). As for the other components, T x , T y , S xx , S yy , and S xy behave like scalar functions, 
and (S xz , S xw ) and (S yz , S yw ) behave like z and w components of vector functions (T z , T w ). 
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C. 4D space with U(l) x U(l) symmetry 



Next, we describe the cartoon method in a 4D space of U(l) x £7(1) symmetry, where 
two Killing vector fields, d v = xd y — yd x and = zd w — wd z , are present. An example 
of a spacetime of this symmetry is a 5D rotating black hole spacetime with two rotation 
parameters [231]. In this spacetime, the black hole is rotating with respect to the (x,y)- 
and (z, w) planes simultaneously. Of course, a 5D rotating black hole with one rotation 
parameter is also the case for this symmetry. We refer to such a case as u x = y, z = w" in 
short. 

The spacetime of this symmetry can be simulated as a 2+1 problem in the cartoon 
method. We prepare grid points on the (x, z) plane and a few neighboring grid points in 
both y and w directions. In this symmetry, two cartoons are necessary: The first cartoon 
to generate the data in the y- direction, and the second cartoon to generate the data in the 
w- direction. The symmetric relations for each cartoon are essentially the same as those in 
the previous two subsections. 

In numerical simulation, we have the data of points (x, 0, z, 0) at each time step. Then, we 
apply the first cartoon to generate the data for grid points (x, y, z, 0). After that, we apply 
the second cartoon to generate the data for grid points (x,y,z,w). Then, all the necessary 
derivatives with respect to y and w can be taken and the data can be evolved to the next 
time step. This method may be called the double cartoon method. As we demonstrate in 
Sec. IV, the double cartoon method works well as the single cartoon method. 

D. 4D space with SO (3) symmetry 

Finally, we consider a space of a different type of symmetry, SO (3) symmetry, in which 
the three Killing vectors, $, 1 := yd z — zd y , £ 2 := zd x — xd z , and £ 3 := xd y — yd x are present. 
In other words, each hypersurface of w = const, is spherically symmetric. An example for a 
spacetime of this symmetry is a head-on collision of two black holes moving along the w axis, 
because the other directions x, y, and z are equivalent if the black holes are not rotating. 
We refer to this case as x = y = z, w in short. 

This spacetime can be simulated as a 2+1 problem in the cartoon method. We prepare 
grid points in the (x, w) plane and a few neighboring grid points in both y and z directions. 
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Using the data on the (x, w) plane, the data at points (x, y, z, w) can be calculated by the 



SO (3) symmetric relations using the data at the point (r, 0, 0, w) where r = ^ x 2 + y 2 + z 2 . 
The data at the point (r, 0, 0, w) should be determined by an interpolation as before. 

In this case, the symmetry relations are different from those in the previous cases. For 
scalar functions, it is trivial as 



ty(x, y, z, w) = \I/(r, 0, 0, w). 



(28) 



In order to derive the symmetry relations for vector and symmetric tensor functions, we 
have to know the S'0(3)-symmetric forms of a vector and a symmetric tensor, which can be 
found by the conditions C^ n T l = and C^Sij = 0. For this purpose, we first describe their 
components in the spherical-polar coordinates (r, 9, <p, w) introduced by re = r sin 9 cos <p, 
y = r sin 9 sin ip, and z = rcos9. Then, in the SO (3) symmetry, they are written as 



T(r,9,if,w) = (T r (r,w), 0, 0, T w (r,w)), 



(29) 



and 



S i:i (r,6,<p,w) 



(30) 



^ S rr (r,w) S rw (r,w)^ 

* S ee (r,w) 

* * Sge(r, w) sin 2 9 
V * * * S ww (r,w) J 

Now we transform these expressions to the Cartesian coordinates and use the relations on 
the (x, w) plane to give 

T x (x,y,z,w) = (x/r)T x (r,0,0,w), (31) 
Ty(x,y,z,w) = (y/r)T*(r,0,0,w), (32) 
T z (x,y,z,w) = (z/r)T x (r,0,0,w), (33) 

for a vector function and 



S xx (x,y,z,w) = (x/r)S xx (r,0,0,w) + (l-x/r)Syy(r,0,0,w), 
S yy (x,y,z,w) = (y 2 /r 2 )S xx (r,0,0,w) + (1 -y 2 /r 2 )S yy (r,0,0,w), 
S zz {x,y,z,w) = (z 2 /r 2 )S xx (r,0,0,w) + (1 - z 2 /r 2 )S yy (r, 0,0, w), 

S yz (x,y,z,w) = (yz/r 2 )[S xx - S yy ](r,0,0,w), 
S zx (x,y,z,w) = {zx/r 2 )[S xx - S yy \{r,0,0,w), 
S xy (x,y,z,w) = {xy/r 2 )[S xx - S yy \{r,0,0,w), 



(34) 
(35) 
(36) 

(37) 
(38) 
(39) 
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for a tensor function. Note that T w and S ww satisfy the symmetry relation of a scalar 
function, and (S wx , S wy , S wz ) satisfy that of (x, y, z) components of a vector function. 

Using these symmetry relations, one can calculate the data at grid points in the neigh- 
borhood of the (x, w) plane (i.e., y,z 7^ 0), and thus the evolution can be performed as a 
2+1 problem. Note that by eliminating the w direction, the above symmetry relations can 
be used also for simulating a 3D spherically symmetric space in a 4D spacetime. 



IV. CODE TESTS 



In the previous two sections, we have described necessary ingredients for higher- 
dimensional numerical relativity, i.e., the BSSN formalism and the cartoon method. Based 
on these, we have implemented several codes for simulating 5D spacetimes in the following 
manner. As often done in the 4D numerical relativity (e.g., Ref. 37|). we adopt the centered 
fourth-order finite differencing in the space directions, except the advection terms such as 
/3 dkjij for which the fourth-order upwind finite differencing is adopted. The time evolution 
is carried out using the fourth-order Runge-Kutta method, where the Courant number is 
adopted to be 0.5. Vertex-centered grids are employed for all the space directions. In the 
present codes, we do not implement adaptive mesh refinement (AMR) algorithm. We plan 
to combine our codes with our AMR code (SACRA code) in the future [381 ] . 

So far, we have developed the 3D codes for spacetimes with U(l) symmetry (x, y, z — w), 
the 2D codes for spacetimes with U(l) x U(l) symmetry (x = y, z — w) and with SO (3) 
symmetry (x = y = z, w). The authors of this paper have independently developed the 
codes, and checked that the numerical results for test simulations derived by the two codes 
agree. In addition, Yoshino has made a ID code for spacetimes with SO (4) symmetry 
( u x = y = z = w" in short). In the following, we present the results by Yoshino's code, 
for which the uniform grid with the grid spacing Ax is always adopted for all the space 
directions. 

In order to prove the validity of our codes, we consider that at least the following two test 
simulations have to be successfully carried out as in the 4D numerical relativity. One is the 
evolution of the 5D Schwarzschild black hole, and the other is the evolution of a spacetime 
composed of gravitational waves of small amplitude. Since their metrics are analytically 
given, they can be used in the benchmark tests for calibrating the codes. The results of 
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the test simulations are reported in Sees. IVA and IVB, respectively. In addition, we show 
in Sec. IVB that the energy flux can be correctly calculated for linear gravitational waves 
using the Landau-Lifshitz pseudotensor. 



A. 5D Schwarzschild spacetime 



First, we analytically derive the 5D Schwarzschild metric in the geodesic slices and then 
compare numerical results with it. Next, we demonstrate that a long-term evolution of the 
5D Schwarzschild spacetime is feasible in the so-called puncture approach, as in the 4D case 



e.g., Refs. |Sl37|,|4l|,|4j). 



1. Geodesic slices 



The well-known metric of a 5D Schwarzschild spacetime is 



ds 2 



-f(r)dt 2 + 



2 1 dr " +r 2 dn 



f(r) 



f(r) 



,2 ' 



where dVt\ is the line element of a 3D unit sphere and is the Schwarzschild radius 



8G 5 M 



3tt 



(40) 



(41) 



Here, we consider the Gaussian normal coordinates starting from the t = hypersurface, 



43 



44]. 



which is analogous to the Novikov coordinates in the 4D Schwarzschild spacetime 
Let us introduce a geodesic congruence of test particles that are initially at rest. Each 
geodesic labels the radial coordinate. Denoting r as the initial value of r for each geodesic, 
we define f by 



r 



f|1 + 4# 



(42) 



as the radial coordinate. At t = 0, the coordinate f agrees with the so-called isotropic 
radial coordinate. Adopting the proper time r for each geodesic as the time coordinate, the 
geodesic equations are solved to give 



2 J2 



r =r - (r h /r ) r 



and 



t = v/TRr + ^log 



(roAh)vT(ro) 



r - ( r i/ r h)Vf( r o) 



(43) 



(44) 
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FIG. 1: Left-hand panel: Snapshots of % x along the x axis for r/r^ = 0.5, 0.6, 0.7, 0.8, and 0.9. 
The unit of x is The grid resolutions are Ax = 0.1 (x) and 0.05 (©). The solid curves denote 

(a) 

the analytic solutions, 7^ '. Right-hand panel: The averaged error, 5 7 , as a function of Ax. Here, 
the average is taken for the data in the range < x < 5. The upper short line segment shows the 
relation of the fourth-order convergence (i.e., a segment with the slope 4). 

Substituting these equations into Eq. (flUI) and transforming r to f with Eq. (j4"2l . we obtain 
ds 2 = -dr> + ^^Hf ^ + H - (W^of ^] ^3- (45) 

bo - (W r o) r ] h 

This line element shows that the ff component of the metric diverges at r = r^/r^. Curva- 
ture invariants indeed show that the curvature singularity appears at this time. This implies 
that at the time r = r^, the slice hits the singularity at f = 

The derived line element (1451) shows the exact solution, and thus it can be used for test 
simulations. In this test, we perform a simulation with the gauge condition a = 1 and 
f3 l = 0, until the computation crashes approximately at the crash time r crash = r h . 

The left-hand panel of Fig. [T] shows the snapshots of xx component of the conformal 4D 
metric j xx along the x axis for various time slices as r/r^ = 0.5, 0.6, 0.7, 0.8, and 0.9. For 
this plot, the grid resolutions Ax = 0.1 and 0.05 are adopted. Here, the units of x are 
(i.e., the event horizon is initially located at x — 1). We see that the values of j xx rapidly 
increase and blow up around x = 1, and agree approximately with the analytic solutions 
fl4"5]) (solid curves). 

The right-hand panel of Fig. [1] plots the averaged error as a function of the grid spacing 
Ax. Here, the averaged error is defined by 

5^ = ^jy x \% x -tll (46) 
17 




FIG. 2: Left-hand panel: The violation of the Hamiltonian constraint at the time r/r^ = 0.9 for 
the grid resolutions Ax = 0.1 (x), 0.05 (0), and 0.025 (•). Although the violation grows as the 4D 
hypersurface approaches the singularity, it becomes smaller for a fixed value of r as the resolution 
is increased. Right-hand panel: The average, 8h, as a function of the grid spacing Ax. The upper 
short line segment shows the relation of the fourth-order convergence (i.e., a segment with the 
slope 4). 

where ^ xx and 7^ are numerical and analytic solutions, respectively, and the integral is 
numerically performed using the data on the grid points. This figure indicates that the 
numerical error approaches zero approximately at the fourth-order convergence. 

Here, the figures are plotted for the results obtained by the ID code (x = y = z — w), 
but essentially the same results are obtained by the 2D codes (x = y, z = w and x = y = z, 
w) and the 3D code (x, y, z = w). 

The left-hand panel of Fig. [2] plots the violation of the Hamiltonian constraint along the 
x axis at r/r^ = 0.9. Here, the violation is defined by 

H° :=R- ^K 2 + AijA'v. (47) 

As the surface approaches the singularity, the value of the constraint violation grows rapidly. 
However, if we fix the time r and compare the results by the different grid resolutions 
Ax = 0.1 (x), 0.05 (©), and 0.025 (•), the clear convergence is seen. 

The right-hand panel of Fig. [2] plots the averaged constraint violation 5h- Here, the 
average is defined in the same manner as Eq. f|46|) . This figure indicates that the error 
converges also at the fourth order approximately. 
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FIG. 3: The evolution of the trace of the extrinsic curvature K at x = (y = z = w = 0) in 
the puncture gauge. The unit of t is The value of K asymptotes to zero. 

2. Long-term evolution 



Next, we show that long-term evolution of a black hole is feasible using the puncture 
approach as in the 4D case. In this test, the initial condition is prepared in the isotropic 
coordinates, and then, the evolution is carried out without excising the black hole interior. 
As the gauge conditions, we adopt the generalized version of the dynamical slicing condition 



39, 



-r] a aK, 



and the T-driver gauge condition [39J] 

(D-l) 



d t (3 l 



2(D 



—ir R l 



(48) 



(49) 



Here, vi ong indicates the propagation speed of a gauge mode and has to be chosen as < 
viong < 1- We tried the choices vi ong = 1 and Vo/2, and found that the stable numerical 
evolution is possible in both cases. The choice vi ong = v^3/2 stabilizes the numerical evolution 
near the puncture a little more. t] a and r]p are positive constants that can be arbitrarily 
chosen. For rj a , we chose several values between 1.2 and 2.0, and found that the stable and 
long-term simulation is feasible irrespective of the value of r\ a . For r]p, we choose l/5r^. 



In the following, we show the results of the numerical evolution for the case vi 



ong 



1 and 



r] a = 2. The initial condition of the lapse and shift is chosen as a = ^fx an d ft = 0. Figure [3] 
shows the evolution of K at x = rhjl on the x axis. The unit of the length is The value 
of K relaxes to zero after a few oscillations, and the slice asymptotes to a maximal surface 
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a 





FIG. 4: The values of the lapse a and the shift vector (3 X along the x axis at the time t = 50r^. 
Here, the unit of x is 



because of the property of the dynamical slicing condition We evolved this spacetime 
up to t = 100r/j, and the spacetime relaxes to a stationary state. Figure H] shows the values 
of a and /3 X along the x axis at t = 50?v By this time, the variables approximately relax 
to stationary states. These results are quite similar to the evolution of a 4D Schwarzschild 
spacetime (compare with Figs. 1 and 5 in Ref. 42j). 

The left-hand panel of Fig. [5] plots the violation from the Hamiltonian constraint defined 
by Eq. (|47p along the x axis at the time t = 50r^ for the grid resolutions Ax = 0.1 (x), 
0.05 (©), and 0.025 (•). After the long-term evolution, the violation in the neighborhood 
of the puncture x = grows large to become 0(1). This is because the analyticity of the 
solution is violated at the puncture. However, the error rapidly decreases as x is increased, 
indicating the reliability of the numerical results. It is also found that the spatial patterns 
of H° depend on the resolution after the longterm evolution, t ^> r^, although initially they 
have similar shapes. 

The right-hand panel of Fig. plots the averaged constraint violation 5h in the range 
0.5 < x < 10 defined in the same manner as Eq. ()46|) . Because of the error generated at 
the puncture, the value of Sh does not show the fourth-order convergence. Nevertheless, the 
violation rapidly decreases as the grid resolution is increased. 

The obtained stationary data are expected to agree with those of the limit surface of 
the maximally sliced evolution. In the 4D case, the limit surface of K = was analytically 
determined [45|, |46( and also the asymptotic solution in the numerical simulation agrees 



with it. A simulation 



421 ] also demonstrates that the spacetime remains in a stationary 



20 



-5.5 



2 4 6 8 10 _i,7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 

x log 10 Ax 

FIG. 5: Left-hand panel: The violation of the Hamiltonian constraint at the time t = 50r/j for the 
grid resolutions Ax = 0.1 (x), 0.05 (0), and 0.025 (•). After the long-term evolution, the spatial 
pattern of H° depends on the resolution, and the error generated at the puncture is 0(1). But the 
general tendency is that the violation becomes smaller as the resolution is increased. Right-hand 
panel: The average, 5h, in the range 0.5 < x < 10 as a function of the grid spacing Ax. The lower 
short line segment shows the relation of the fourth-order convergence (i.e., a segment with the slope 
4). The convergence is worse than the fourth-order convergence because of the error generated at 
the puncture. 

state if the limit surface is adopted as the initial condition. The limit surface exists also 
in a higher-dimensional Schwarzschild spacetime and it provides a useful benchmark for 
calibrating codes for higher dimensions, as shown in Ref. [471 ]. 



B. Linear gravitational waves 

We turn our attention to a simulation for propagation of gravitational waves of small 
amplitude. Here, we focus only on gravitational waves that preserve U(l) x U(l) symmetry 
(x = y, z = w) or 5*0(3) symmetry (x — y — z, w). In Appendix A, the linearized Einstein 
equations of such symmetries and their special solutions for the lowest multipole moment I 
are described. In this subsection, we pick up a tensor-mode perturbation with U(l) x U(l) 
symmetry for the test simulation. 

As the perturbative solution used for the test simulation, we adopt the spatial metric 
(IA35j) with Eq. (IA34I) and the special solution for h(t, r) given by Eq. (IA23I) with the gauge 



condition a = 1 and = 0. Here, we set A = 1/6 and B = in Eq. (TA"3~4"]) . and A = 0.015 
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FIG. 6: Left-hand panel: Snapshots of ^ zz along the x axis for t = 1,2, 10 for propagation of 
a gravitational-wave packet. The dotted points (•) and solid curves denote the numerical results 
and perturbative solutions, respectively. Right-hand panel: The averaged deviation of j xx from 
the analytic perturbative solution (■) and that from the numerical data with the grid resolution 
Ax = 0.05 (□) as functions of Ax. Here, the data at time t = 3 are used and the average is taken 
for the data in the range < x < 6 and < z < 6. The upper short line segment shows the 
relation of the fourth-order convergence (i.e., a segment with the slope 4). 

and uq = 2 in Eq. (1A23I) . In the simulation, we evolved the initial data that correspond to 
the perturbative solution under the same gauge condition a = 1 and (3 l = 0. The left-hand 
panel of Fig. [6] compares the analytic solution for the linearized Einstein equation (solid 
curves) and the numerical results (dotted points) obtained by the 2D u x = y, z = w" code 
(where the double cartoon method is used). The values of j zz are plotted along the x axis 
for t = 1,2,..., 10. The two results agree well, indicating the validity of our code. 

The right-hand panel of Fig. [6] shows the deviation of the numerical solution from the 
analytic solution of the perturbation as a function of grid sizes Ax (black squares, ■). Here, 
we used the data of ^ xx on the (x, z) plane at t — 3 and evaluated the deviation by taking the 
average of \"y xx — 7^ | in the region < x < 6 and < z < 6. The deviation scarcely depend 
on Ax, and thus it is not caused by the grid resolution. The deviation primarily comes from 
the fact that the perturbative solution ignores the second- and higher-order quantities in hij, 
whereas the numerical simulation is carried out by the fully nonlinear evolution equations. 
Indeed, the order of the difference ~ 10~ 5 agrees with the magnitude of the nonlinear effect 
for our chosen wave amplitude. 

Squares (□) in the same panel show the difference of the numerical data computed with 
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the grid resolutions Ax > 0.1 from the one with the grid resolution Ax = 0.05. The 
difference decreases approximately at the fourth order, implying that our numerical solution 
achieves the fourth-order convergence. 

Using the analytic solution of the linearized Einstein equation, we can also test extraction 
methods of gravitational-wave energy flux from the numerical data. The Newman-Penrose 
variable is now widely used for extracting radiated energy of gravitational waves in the 4D 
numerical relativity. However, the formalism for the extraction based on the Weyl scalar 
has not been developed in higher dimensions. Thus, for the calculation of the energy flux, 
we adopt the Landau-Lifshitz pseudotensor C U whose D- dimensional form is given m 

LL 



Appendix B. Because is not a tensor, it is a coordinate-dependent quantity in general. 



However, as discussed in p. 85 of Ref. 



19] , the total amount of gravitational energy and the 



total radiated energy, Eqs. (1B8I) and (1B9I) . are shown to be the gauge- invariant quantities for 
the linear gravitational waves of a perturbed flat spacetime up to second order with respect 
to the metric perturbation. 

Let r = r Q bs be the radius of the extraction of gravitational-wave energy flux. The 
integrated energy flux E TSbd through the surface r = r obs is evaluated by 

£ rad (r obs ) = J t&dSdt, (50) 

where dS is the area element of a hypersphere of r = r Q ^ s . Here, the second-order expression 
of the Landau-Lifshitz pseudotensor with respect to the perturbative quantities is used (see 
Eq. flB6j) and subsequent explanation). We evaluate -E ra d(^obs) both for the perturbative 
solution and for the full numerical solution, and compare the two results. 

The value of -E" r ad( r obs) for the perturbative solution is evaluated semianalytically by 
proceeding the integrations of Eq. (150)) numerically. Figure [7] shows the value of StiG^E^uIq 
as a function of r Q b s ^o by the solid curve. E rad changes from zero to a constant value 
~ 0.00925/87rG5u;o as r Q ^ s uo increases from zero to 40. The value of E iad near the center 
^obs^o ~ 1 does not have a definite meaning because the Landau-Lifshitz pseudotensor is 
not gauge invariant. However, the asymptotic value of E rad for r obs u; 3> 1 is gauge invariant 
and should indicate the correct amount of the radiated energy. Note that because of the 
conservation law <\B7\i . d^t^[ = 0, the integrated energy flux -E ra d( r obs) has to be equal to 
the initial amount of energy within the surface r = r Q b s , 

E(r ohs ) = I t™ L dV (51) 
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FIG. 7: The integrated energy flux SttG^E^Uq calculated by the Landau-Lifshitz pseudotensor as 
a function of the extraction radii r b s wo for the analytic solution (solid curve) and for the numerical 
results computed with the grid resolutions Ax = 0.1 (0), 0.15 (a), and 0.2 (x). The inset shows 
the enlarged figure for the region 31 < r b s wo < 40. 

where dV is a volume element of the 4D space. This is directly checked by the numerical 
integrations. 

In order to calculate -E'radtjobs) for the numerical solution, we proceed as follows: At each 
time step, we define the perturbed quantities as h u = — 2(a — 1), h ti = fy, and hy = 7ij — <%, 
and evaluate t l ° at r = r b s using Eq. (IB6I) . Then, we calculate the integral 



that gives the energy flux through a surface r = r b s , where hi is the outward unit normal to 
the surface. Finally, dE Ta ^/dt is integrated from the initial time to the final time to obtain 
-E'rad(^obs)- In order to evaluate the metric functions for a surface of r = r b s , we used linear, 
quadratic and cubic interpolations and compared the results. Although relatively large 
deviation from the analytic result is seen for the linear interpolation, the deviation becomes 
smaller when the quadratic interpolation is used. The result of the cubic interpolation did 
not improve that of the quadratic interpolation. This is because in these cases, the error 
primarily comes from the error generated at the outer boundary (i.e., the error due to the 
inaccuracy of the outgoing boundary condition) which is approximately at the second order 
with respect to the grid size. 




(52) 
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Figure [7] shows the results of SttG^E^^o evaluated by the quadratic interpolation at 
several observation points for the grid resolutions Ax = 0.1 (©), 0.15 (A), and 0.2 (x). 
The deviation decreases as the resolution is increased and the numerical data approach the 
analytic result. 

We also evolved the gravitational- wave packet using the dynamical slicing and T-driver 
conditions and checked that _E ra d depends very weakly on the initial choice of the lapse and 
shift as long as the initial value of a — 1 and f} 1 is small. This is natural because the gauge 
invariance of E rad is guaranteed for r obs c<; 1. Therefore, we conclude that the extraction 
of gravitational-wave energy flux by the Landau-Lifshitz pseudotensor works well, as far as 
the amplitude of gravitational waves is sufficiently small at the extracted region (i.e., in the 
wave zone). 

V. SUMMARY 

This paper describes the formulations for numerical relativity in higher dimensions and 
reports the new codes for simulating 5D spacetimes. We derived the BSSN formalism for 
higher-dimensional spacetimes and also studied the cartoon method in 5D spacetimes of U(l) 
symmetry (x = y,z,w), U(l) x £7(1) symmetry (x — y, z — w), and 5*0(3) symmetry (x = 
y = z, w). Based on the BSSN formalism and the cartoon methods, we have implemented 
the new 5D numerical relativity codes, and tested them by evolving the 5D Schwarzschild 
spacetime and the spacetime composed of gravitational waves of small amplitude. The 
numerical results converge to the analytic solutions with improving the grid resolution at the 
correct order (fourth order). It was also demonstrated that the 5D Schwarzschild spacetime 
can be evolved for a long time by the puncture approach, as in the 4D case. 

We proposed the method of extracting gravitational-wave energy flux by the Landau- 
Lifshitz pseudotensor and tested this method. We showed that the integrated energy fluxes 
calculated at several surfaces r = r b s agree well with the semianalytic solution derived by 
perturbative calculations. Furthermore, it was confirmed that the result is insensitive to the 
gauge conditions for the lapse and shift. These results indicate that the energy extraction by 
the Landau-Lifshitz pseudotensor works well. The remaining problem would be to check that 
the extraction in the wave zone works well even when the central region is highly nonlinear, 
e.g., a head-on collision of two black holes. This will be tested by performing a simulation of 
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Brill wave spacetime in five dimensions and by comparing the ADM mass of the initial data 
and the energy radiated during the evolution. Also, it is necessary to check if the extraction 
of the angular momentum is possible. We expect that the radiated angular momentum also 
can be calculated by the Landau-Lifshitz pseudotensor using a similar manner to the 4D 
case [481 ] . 

As discussed in Sec. I, there are many interesting issues of nonlinear dynamics in higher- 
dimensional gravity, which should be studied in numerical relativity. In this paper, we have 
prepared the tools necessary for simulating higher- dimensional spacetimes. Our next step is 
to apply our codes to the unsolved problems. 
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APPENDIX A: LINEAR GRAVITATIONAL WAVES IN THE 5D SPACETIME 



In this section, we describe solutions of the linearized Einstein equations in the 5D flat 
spacetime focusing on the perturbations preserving £7(1) x U(l) symmetry and SO (3) sym- 
metry. In the following, we denote the metric perturbation as h^, which obeys the linearized 
Einstein equation 

5GV[M = 0. (Al) 
The analysis for perturbations of higher-dimensional Schwarzschild black holes described 



in Ref . 



491 ] is partially used, since this formulation is applicable also for the flat spacetime. 



In their approach, the perturbation is decomposed into the scalar, vector, and tensor modes 
(with respect to the 3D unit sphere) using spherical harmonic functions, and the master 
equations are derived for the gauge-invariant variables. We adopt their method for spherical 
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harmonic expansion but do not use the master equations, because we are interested in 
explicit special solutions for which the master equations are not necessary. 



1. Perturbation with U(l) x U(l) symmetry 

We here derive a solution of U(l) x U(l) symmetry for Eq. flAlj) . Denoting the coordinates 
by (t,x,y, z,w), we introduce the following curvilinear coordinates: 

x = r sin 9 cosip, (A2) 

y = r sin 9 sin (p, (A3) 

z = r cos 9 cosip, (A4) 

w = r cos 9 sinip. (A5) 

In these coordinates, the line element of the flat space is 

dl 2 = dr 2 + r 2 (d9 2 + sin 2 9dy 2 + cos 2 9di\) 2 ). (A6) 

For the U(l) x U(l) symmetry case, the linear perturbation h^ v satisfies 

dh dh 

-t = W = ' (A7) 

since d v and are the Killing vectors. 

In the gauge condition with a = 1 and [3 k = (i.e., h 00 = h 0i = 0), the spatial components 
of the linear perturbation hy satisfy 

kj = &h ij} (A8) 

where A is the flat 4D Laplacian. The Hamiltonian and momentum constraints in the linear 
approximation give 

f%5 = 0, (A9) 
D%j = aj , (A10) 

where is the 4D flat space metric in the curvilinear coordinates and Di is the covariant 
derivative with respect to g^. aj denotes a constant vector determined at the initial state, 
which is set to be zero in the following for simplicity. 
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a. Scalar mode 



The scalar mode (with respect to a 3D unit sphere) is expanded in terms of scalar har- 
monic functions § on a 3D unit sphere with the metric da 2 = ajjdz 1 dz J = d9 2 + sin 2 9d<p 2 + 
cos 2 9dtp 2 , which satisfy the equation 



[A3 + 1(1 + 2)]S = 0, 



(All) 



where A 3 = D 1 Dj is the Laplace operator on the 3D unit sphere. In the following, we 
focus only on solutions of the lowest-order multipole moment I = 2, for which the harmonic 
function is 

S = 2cos 2 0-1. (A12) 



The scalar-mode perturbation is given in the form 



a(t,r) 



rb(t, r) 



where 



§j := Djt 



* ^[c^r^au + dfar^u] 



and ^u^j^DtDjS + lvul 



or more explicitly, 



i/J = - x 




\ 



§j = ( — 4 sin 9 cos 9, 0, 0), 

1 - 2 cos 2 9 

* sin 2 fl(cos 2 fl-2) 

* * cos 2 9(1 + cos 2 9) J 



Equation flA~9j) gives a + 3c = 0, and Eq. flAlOl) yields 



8b = 4a + ra , 



—d = 46 + r& r + c. 



The rr component of Eq. flA8j) with Eq. (IA17I) gives a wave equation for a: 

7 

Ql Ct yr> | CL f m 



(A13) 

(A14) 

(A15) 
(A16) 



(A17) 
(A18) 

(A19) 



Equations (IA17I) and (IA18[) imply that once a is computed from the equation (IA19I) . 6, c 
and d are subsequently determined. The obtained solution is guaranteed to satisfy other 
components of Eq. (1A8j) . 
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Defining a = u/r 3 , we obtain the equation 



1 n 2 

U = U trr H U t r —It. (A20) 



Here, n — I + 1 = 3. The formal solution of this equation is written as 



u = Re 



dujf{uj)e l ^J n {uor) 



(A21) 



where f(co) is an arbitrary function of lu, and J n is the Bessel function of n-th order. In the 
integral expression, it is written by 

i r 2n 

Jj z ) = — / dticos{n$ - zsintf). (A22) 
27r Jo 

To constitute a solution for the propagation of a gravitational-wave packet, we set f(u) = 
-i^/2nA e~ w2/2ul o. Then, Eq. (IA2ip is integrated to give 

u(t,r)=A uJo ^sin(n^)e-^ ( *- rsin,5)2/2 . (A23) 
Jo 

In this case, u = at t = 0, and thus fey = 0, whereas the extrinsic curvature Kij = —hij/2 
is not zero because 

-2- 



U 



(0,r) = rA uj 3 dtism(n<&) sm^e- {uJorsin ^ )2/2 (A24) 
Jo 



is not zero at t = 0. 



6. Vector mode 



Perturbation foy of the vector type is expanded in terms of the harmonic vectors V/ 
satisfying 

[A 3 + /(Z + 2) - l]Vj = 0, (A25) 

DjY J = 0. (A26) 

Under the U(l) x C/(l) symmetry, only the modes for odd I numbers are nonzero. Since the 
/ = 1 mode denotes a stationary perturbation with angular momentum, the lowest value of 
I is 3 for the nonstationary perturbation. For this mode, 

V J = (0, A(sin 2 6 - 2/3), ^(cos 2 9 - 2/3)) , (A27) 
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where A and B are arbitrary constants. The perturbation is given in the form 



Kl 



(l/r)k(t,r)Vj 
* ro(t,r)Y u 



where V/j is defined by 



/ 



V 



DjYj + DjVj 



,4 sin 3 6» cos 6» -B sin 6 cos 3 6 » 








From Eqs. (IA8p and flAlOp . the equations for /c and o are derived as 

v_ 1 16 



6o = /c r H — A;. 

r 



(A28) 



(A29) 



(A30) 



(A31) 



Here, Eq. flA30|) has the same form as Eq. (1A20|) but with n = 4. Hence, a special solution 
for k(t,r) is given by Eq. (I A23[) with n = 4, and then o(t,r) is calculated from Eq. (IA31I) . 



c. Tensor mode 

Perturbation h%j of the tensor type is expanded in terms of the harmonic tensors Tjj 
satisfying 

[A 3 + 1(1 + 2) - 2] Tjj = 0, (A32) 
T*j = 0, DjT J j = 0. (A33) 
Under U(l) x £7(1) symmetry, the possible harmonic tensors for / = 2 are 

( A ^ 

T 7J = * A sin 2 0(1 -3 sin 2 6) B sin 2 6 cos 2 9 

^ * * A cos 2 8{3 sin 2 - 2) y 

where A and 5 are arbitrary constants. The perturbation is given in the form 



hij 





* rh(t,r)Tu 



(A34) 



(A35) 
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and from Eq. (1A8I) . the equation for h becomes 



19 

h = h^ rr H — h :r -h. (A36) 



This is the same equation as Eq. (IA20[) . and thus a special solution is given by Eq. (1A23P 
with n = 3. 

2. Perturbation with SO(3) symmetry 

Next, we derive solutions of a perturbation of 5*0(3) symmetry. For the Minkowski 
coordinates (t,x,y,z,w), we introduce the following hyperspherical coordinates: 

x = r sin 8 sin (p sin ip, (A37) 

y = r sin^sin^cos^, (A38) 

z = r sin 9 cos (p, (A39) 

w = rcos6. (A40) 

Then, the line element of the flat space is given by 

dl 2 = dr 2 + r 2 (d6 2 + sin 2 6dcp 2 + sin 2 9 sin 2 cpdip 2 ). (A41) 

Here, we consider solutions of 50(3) symmetry with the Killing vectors 

£i = — cos ifjd^ + cot if sin ipd^, (A42) 

£ 2 = sinipdip + cot ip cosipd^, (A43) 

^ 3 = ^. (A44) 

Under the requirement of this symmetry, the vector and tensor modes do not exist, because 
there are no vector and tensor harmonic functions that satisfy £f n Vj = and C^ n Tjj = 0. 
Therefore, only the scalar mode should be analyzed. 

The scalar harmonic function defined by Eq. (IA11I) on a 3D unit sphere with the metric 
da 2 = a IJ dz 1 dz J = d9 2 + sin 2 9dip 2 + sin 2 9 sin 2 Lpdip 2 is 

§ = 4 cos 2 9-1, for / = 2. (A45) 



The metric perturbation is given in the same form as Eq. (1A13|) . Here, definitions for §/ 



and §/j are same as Eq. flA14p . and their explicit forms in this case are 

§ J = (-8 sin 9 cos 9, 0, 0), (A46) 
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hj = - sin 2 6 x 



\ 



V 



(A47) 



2 

* -sin 2 # 

* * — sin 2 6 sin 2 ip 
The equations for a, b, c, and d are the same as Eqs. (1A17I) - (1A19|) . and a special solution 
for u = r 3 a is given by the same formula as Eq. (IA23I) . 



/ 



APPENDIX B: LANDAU-LIFSHITZ PSEUDOTENSOR 



In this section, we derive the Landau-Lifshitz pseudotensor in a D- dimensional spacetime 
M. with the metric g^. Following [40J, we define 



-99 



(Bl) 



where g is the determinant of the metric, and then we introduce the super-potential 

H l**V0 = g^g*P _ g^grf. (B2) 



The Landau-Lifshitz pseudotensor is defined by 

WirGntZ = (- g )-^H^ ap - (2itT - <TR) . 
From this definition, the conservation law is derived: 



(B3) 



(B4) 



Because the Landau-Lifshitz pseudotensor is not a tensor, it does not have a covariant 



meaning in general. However, for a perturbed flat spacetime, the leading-order terms of t 



fiV 

LL 



with respect to the perturbative quantities can be used to evaluate the total energy and 
total radiated energy of the gravitational field in a gauge- invariant manner (see below). 



In Ref. 



40j, two expressions for t^ h are given. The first one is the expression in terms 



of the Christoffel symbols, and this expression holds for arbitrary dimensionality D. The 
second one is the expression by the metric functions, and it is modified to give 

i6nG D (-g)tZ = r:j aP p - r a ,J v % + \g» u g aP ~g ap J^ P 

- {g^g Pp ? p ja f° a + g ua g Pp r p ,J^ a ) + gapsTsTjr^ 
l 



+ 



(2</"V - gTsf*) [(D - 2)g pa g jS - g a ,g p5 ] g^J*^ (B5) 
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in D dimensions. Let us consider the perturbation on a flat spacetime, whose metric is 
9y,v = Vfj,f + hfxv, where T]^ is the flat metric in the Minkowski coordinates. Defining : = 
— (l/2)hrj liu , the Landau-Lifshitz pseudotensor is rewritten as 



a,p 

+ \h^h p f - ^h^ a h pa>a - (2W - r u h' a h a } . (B6) 

Here, we have kept only the second-order quantities of the perturbation. Note that the 
second-order Landau-Lifshitz pseudotensor behaves as a tensor against the general co- 
ordinate transformations of the background spacetime (but for a fixed gauge), by replacing 
the coordinate derivatives to the covariant derivatives and the Minkowski metric 7] pl/ 
to the flat background metric g pv in curved coordinates in Eq. (IB6I) . The quantity t 0r in 
Eq. (|50|) has to be evaluated in this way. 

For the expression flB6p . the conservation law flB4[) for a vacuum spacetime becomes 

<Vr L = (B7) 

in the Minkowski coordinates, which suggests that i^L can be interpreted as the effective 
stress-energy tensor of the gravitational field valid up to second order in hn U . Here, we have 



to be careful because the Landau-Lifshitz pseudotensor is not the unique quantity satis 
the condition (1B7|) and also because this quantity is not gauge invariant (see Ref. 
However, the total energy 



ymg 



id) 



£tot = J t™ L dV (B8) 

is shown to be the gauge-invariant quantity, where dV is the volume element of the spacelike 
hypersurface. Similarly, the total radiated energy 

£ rad = J t 0i hidSdt (B9) 

is gauge-invariant, where dS and n l are the area element and an outward unit normal of a 
surface at the distant region. Therefore, the Landau-Lifshitz pseudotensor ££l provides us 
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a reliable method for evaluating the total radiated energy. 
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